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Abstract 

The spontaneous self-organization of two-dimensional magnetized plasma is investigated within 
the framework of magnetohydrodynamics with a particular emphasis on the symmetry-breaking 
induced by the shape of the confining boundaries. This symmetry-breaking is quantified by the an- 
gular momentum, which is shown to be generated rapidly and spontaneously from initial conditions 
free from angular momentum as soon as the geometry lacks axi-symmetry. This effect is illustrated 
by considering circular, square and elliptical boundaries. It is shown that the generation of angular 
momentum in non-axisymmetric geometries can be enhanced by increasing the magnetic pressure. 
The effect becomes stronger at higher Reynolds numbers. The generation of magnetic angular 
momentum (or angular field), previously observed at low Reynolds numbers, becomes weaker at 
larger Reynolds numbers. 

PACS numbers: 52.30.Cv, 47.65.-d, 52.65.Kj 
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I. INTRODUCTION 



Understanding the coupling of a magnetic field with the motion of plasmas or conducting 
fluids is a challenging issue both from a fundamental and an applied perspective. In partic- 
ular the self-organization of the velocity and magnetic fields at large scales is an intriguing 
phenomenon. One example is the dynamo-problem, studying the formation of a large scale 
magnetic field induced and amplified by fluid motion (see for example reference^ for recent 
experimental progress). Another example is large-scale spontaneous toroidal and poloidal 
rotation observed in fusion plasmas, an effect that is beneficial for confinement as it may 
suppress turbulence and radially extended structures. This effect may be related to the 
transition to an improved confinement stated. The absence of this transition might jeopar- 
dize the success of the ITER- project. The understanding of large-scale self-organization is 
therefore a key issue in different branches of physics and deserves detailed investigation. 

An academic example of self-organization is the spontaneous generation of angular mo- 
mentum in two-dimensional hydrodynamic turbulence. This phenomenon was discovered by 
Clercx et al. ^ by considering flow in a square domain. We note that this effect was also 
present, but not recognized as such, in calculations by Pointin and LundgrenA. In circular 
domains it was observed to be absent^*^. In^ it was shown that the strength of the spin-up 
can be controlled by increasing the eccentricity of an elliptic domain. For recent reviews 
on the dynamics of two-dimensional turbulence bounded by walls we refer to^*^ and for an 
explanation of spin-up in terms of statistical mechanics to^*^. 

In a recent work^^, it was shown that this effect is enhanced in magnetohydrodynamics. 
The shape of the boundary which contains a plasma may thus be very important in de- 
termining the dynamics of close to two-dimensional plasma flow. In three dimensions, the 
importance of the shape of the plasma container is far from trivial. Indeed, while in infinite 
cylinders plasma can be retained in a static, quiescent state by the Lorentz force, toroidal 
geometries are shown to induce non-zero velocities due to visco-resistive effects^"—. These 
studies concentrated on steady states in axi-symmetric geometry which could be qualified as 
two-and-a-half dimensional. It is reasonable to expect that the same statement will be true 
in fully three-dimensional non-stationary MHD. That case will be studied in future work. 
Here we will consider the unsteady case, but in two space dimensions. 

In the present work we will extend the investigation presented in^. Wall bounded two- 
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dimensional MHD turbulence will be studied, in which the solid boundaries are taken into 
account by the penalization method^. This method is relatively young and has been applied 
to MHD turbulence only recently^^, so that the present paper, in addition to its physical 
relevance, also constitutes a check of the capability of the method to model the influence 
of walls on high Reynolds number MHD turbulence. We consider simulations in which the 
Reynolds number is increased by approximately two orders of magnitude with respect to the 
previous works^^ii^. We consider three differently shaped confining domains. In addition to 
the square and circular geometry considered in the previous study we consider an ellipse. 
The choice of this geometry is inspired by the work of Keetels et al. ^ and this geometry 
has the particularity with respect to the other two to be non-circular, without the presence 
of sharp corners. The initial conditions are completely free from angular momentum, where 
ir)M a small but non-zero initial angular momentum existed. It is shown that the tendency 
to generate angular momentum becomes stronger at higher Reynolds number in the non- 
axisymmetric geometries, while it is absent in the circular container. Furthermore, the 
tendency to generate angular fields vanishes in the limit of large Reynolds numbers. An 
explanation is given for the vanishing of this magnetic angular momentum. 

The remainder of the paper is organized as follows. In Section II the mathematical model, 
the governing equations and their numerical discretization are described. Numerical results 
are presented in Section HI and finally conclusions and perspectives for future work are 
given in Section IV. 

II. MATHEMATICAL MODEL OF BOUNDED MHD TURBULENCE 
A. Governing equations and boundary conditions 

Direct numerical simulation of high Reynolds number MHD turbulence constitutes a chal- 
lenge for computational physics due to the presence of a multitude of nonlinearly interacting 
spatial and temporal scales. Presently the most efficient method to solve homogeneous tur- 
bulence (both hydrodynamic and MHD) is by pseudo-spectral methods, using fast Fourier 
transforma^i^. The additional complexity induced by the presence of solid walls requires 
advanced numerical methods. Pure spectral simulations have been proposed and applied to 
study wall bounded MHD^i, but their prohibitive complexity for increasing Reynolds num- 
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bers limits their application to flows with a relatively limited range of interacting degrees of 
freedom. 

An efficient method to compute flows in the presence of solid obstacles and walls is the 
volume penalization approach which was introduced by Angot et al— for the Navier-Stokes 
equations and applied to hydrodynamic turbulence in^i^I. This method was extended to 
MHD turbulence in a recent work—. Using this method, efficient pseudo-spectral solvers 
can be used to compute flows which contain solid walls and obstacles, which may even move 
in time^^. 

The governing equations are 



with u the velocity, B the magnetic fleld, p the pressure and j = \/ x B the current density. 
Here u and 1] are respectively the kinematic viscosity and the magnetic diffusivity. The 
last term in the evolution equations for u and B is the penalization term which allows to 
impose the solid boundary conditions. Thus both the fluid-domain and the conflning walls 
are embedded in a 27r-periodic square domain. We consider circular, square and elliptic 
domains. For further details we refer to^. 

The quantities no and Bq correspond to the values imposed in the solid part of the nu- 
merical domain. Here we choose Uq = and Bq = B\\ . Here B\\ is the tangential component 
of B at the wall which is not being flxed at a constant value but being re-computed at each 
time-step. Thus the normal component of the magnetic fleld vanishes at the wall, while the 
tangential component can freely evolve. This conflguration corresponds to an electrically 
conducting fluid or plasma in a container with perfectly conducting walls, coated on the 
inside with a thin insulating layer—. In addition to the normal component of the magnetic 
fleld, the current density can not penetrate into the walls, a property which is automatically 
satisfled for two-dimensional flow since the current density only has a component perpen- 
dicular to the plane of the flow. The mask function x is equal to inside the fluid domain 
(where the penalization terms thereby disappear) and equal to 1 inside the part of the do- 
main which is considered to be a solid. The physical idea is to model the solid part as a 
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porous medium whose permeability e tends to zero^^i^^. For e — )■ 0, where the obstacle is 
present, the velocity u tends to Uq and the magnetic field B tends to Bq. The nature of 
the boundary condition for the velocity is thus no-slip at the wall. 

B. Numerical method 

In the case of two-dimensional fiow (here in the x — y plane) it is convenient to take 
the curl of eq. ([H [2]) to obtain after simplification equations for the vorticity and the 
current density, which become scalar valued (in the 2;-direction) and are perpendicular to 
the velocity and the magnetic field, respectively. The vorticity is defined by uBz = V x n 
and = V X B denotes the current density. Furthermore we define the vector potential 
a = ae^ as S = V x a and the stream function as u = V'^ip = {—dip/dy^dip/dx). We 
discretize the evolution equations of vorticity and current density. 



using a classical Fourier pseudo-spectral method. Terms containing products and the penal- 
ization terms, are evaluated by the pseudospectral technique using collocation in physical 
space. To avoid aliasing errors, i.e. the production of small scales due to the nonlinear 
terms which are not resolved on the grid, we de-aliase at each time step, by truncating 
the Fourier coefficients of u and j using the 2/3 rule. For time integration we use a semi- 
implicit scheme of second order, an Euler-Backwards scheme for the linear viscous term and 
an Adams-Bashforth scheme for the nonlinear terms, see e.g.—. 

C. Initial conditions 

The main goal of the present work is the investigation of the formation of large scale struc- 
tures containing significant angular momentum. We therefore want our initial conditions to 
respect two criteria. In the first place we want them to be free from angular momentum. 



— + u- Vuj = B - Vj + uV^cu 
dt 

--(V X Ixiu-uo)]) Bz 
^ + V\[uxB]-ez) = r]V'j 
-- (V X [x{B - Bo)]) ■ 



(6) 



(5) 
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in the second place we want them to be free from coherent structures. One way to gen- 
erate a zero-angular momentum initial condition is, as described in2^, to take an ensemble 
of a large number of Gaussian vortices equally spaced. Half of the vortices have positive 
circulation and the other vortices have negative circulation. The disadvantage is that the 
initial condition hereby contains coherent structures. A straightforward way to generate an 
initial condition without coherent structures, is to start with Gaussian random noise. The 
absence of phase correlations ensures that no structures are present. We therefore initialize 
both vorticity and current density fields with Gaussian random noise as in^^. The Fourier 
transforms cD and j, where uj{k) = J u{x)e~'^^'^dx, are initialized with random phases 
and their amplitudes yield isotropic energy spectra of the form: 

{g + {k/ko))^ 

where g = 0.98 and fco = jV^t^- This energy spectrum is peaked at the largest scales and 
follows a power law proportional to k~^ at large wavenumbers. The energy spectra are thus 
the same for the magnetic field and the velocity field. The phases of the Fourier-modes are 
however chosen randomly and independently, so that the initial fields are different. The 
corresponding fields u and B are calculated from u and j using the Biot-Savart law. The 
fields contain vanishing cross-helicity UiBidA, with Q the flow domain. The so-generated 
fields are however, in general, not free from angular momentum. We note that this was the 
case in reference [l3|, in which the initial conditions contained a small amount of angular 
momentum. We want to avoid this in the present study in order to be able to answer to 
the question whether it is possible to generate angular momentum when initially none is 
present. 

Before describing how we achieve the generation of initial conditions free from angular 
momentum, let us recall the definition of angular momentum Lu and angular field Lb, 
respectively, 

Lu = / ■ {r X u) dA = -2 ipdA, 
Jn Jn 



e^-{rxB)dA =2 adA, (7) 

Q Jq 



where r is the position vector with respect to the center of the domain. Note that the 
equalities on the right hand side assume that a and ip vanish at the boundary of the fluid 
domain. The angular field integral in terms of the vector potential a has some significance 



for 'reduced' MHD^-. To obtain initial fields with = = 0, we proceed as follows. 
We generate one set of fields ui, Bi with corresponding angular momenta and and a 
second set U2, B2 with corresponding angular momenta and L\. By linear combination 
of these conditions, 

u = ui- -f W2 B = B,- -1^2, (8) 
we get initial velocity and magnetic fields free from kinetic and angular momentum. 



III. NUMERICAL RESULTS 



We investigate in total 63 computations in a square, circular and elliptic domain, the 
latter with an excentricity equal to 0.6. The mechanical Reynolds number and magnetic 
Reynolds number are defined, respectively, as 

7^.^^ (9) 

7J«^^. (10) 

V 

The Reynolds numbers are based on the initial root mean square velocity U, the domain 
size D and the kinematic viscosity u and resistivity rj. The magnetic Prandtl number u/r) 
is unity in all simulations so that both Reynolds numbers are equal and denoted by TZ. In 
the following we will therefore not distinguish between the two Reynolds numbers. Two 
series of computations denoted by A and B were performed at a resolution of 512^ grid- 
points and at Reynolds numbers of the order 10^ and 10^ respectively, performing 10 runs 
for each geometry for each Reynolds number. The third series, denoted by C was performed 
at resolution A^^ = 1024^, at Reynolds number of order 10^. The time is normalized by 



D/ ^JIEuit = 0), D being the typical lengthscale of the fluid domain, i.e. the sidelength of 
the square, the diameter of the circle and the longest cross-section of the ellipse. Parameters 
of the simulations are listed in table [H 



A. Visualizations 



Visualizations of the vorticity lo, the stream-function ip the current density j and the 
vector-potential a are displayed in Figure 1. The displayed results are typical results for 

7 




FIG. 1: (Color online) Visualizations of (from top to bottom) the vorticity u, the stream-function 
ip, the current density j and the vector potential a for the square, circular and elliptic geometries. 
The three columns correspond to (from left to right) to the time instants t* = 3, 3, 2.7 of series B 
for which L„ (Fig. 2) is maximal. The time is normalized by the initial turn-over time. Note that 
the numerical method used in the present work does not impose a zero value of a and "0 at the 
wall of the fluid domain. Thus a constant value was substracted from a and ip at every point in 
the fluid domain to impose this. 

series B. We will first focus on the behavior in the square geometry. It is observed that both 
the velocity-field and the magnetic field exhibit a tendency to generate large-scale structures. 
The current-density shows that the magnetic field-lines of the two main flow-structures are in 
the opposite direction. This is even clearer in the plot of the vector potential. The magnetic 
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TABLE I: Parameters of the simulations of series A, B and C. SU*: number of spin-up. The initial 
kinetic and magnetic energies are Eu{0) = 0.3 and -E'fe(O) = 0.7 respectively for all simulations. 
The penalization parameter e is chosen 5 • 10"'* for all runs. 
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angular momentum is therefore small, since the contributions of both structures cancel 
each other out. Note that the right hand side of equation ([7]) relates the magnetic angular 
momentum directly to the vector potential. 

In contrast, the velocity field displays significant symmetry-breaking, which is directly 
reflected in the stream-function. Both vortices are turning in the same sense, with a strong 
shearing region in between them. Non-zero angular momentum results. Similar observations 
can be made for the elliptic geometry. In the circular geometry it is more difficult to visually 
evaluate the generation of angular momentum. 

B. The influence of the Reynolds number and geometry 

To quantify the extent to which a large-scale swirling structure dominates the flow, we plot 
in Figure [2] the angular momentum in the three geometries for series A and B corresponding 
to Reynolds numbers of order 10'^ and 10^, respectively. Since not all runs present spin-up 
(a flow is defined to spin-up when the amount of angular momentum is greater than 10% 
of the angular momentum of a solid-body having the same initial kinetic energy), we 
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FIG. 2: (Color online) Influence of the Reynolds number on the spin-up: time-dependence of the 
absolute value of the normalized kinetic angular momentum averaged over 10 simulations of 
series A (7^ ~ 10^) and series B (7^ ~ 10'*) for the square, circular and elliptic geometry, from 
top to bottom. Here and in the following the angular momentum is always normalized by Cu{0) 
(and Cb{0) for the magnetic equivalent) corresponding to the angular momentum of a solid-body 
having the same initial kinetic energy. 
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FIG. 3: (Color online) Comparison of series B, 7^ 10"^ (top) and series C, 7^ 10^ (bottom). 
Time-evolution of the angular momentum (left) and angular field Lb (right) in the square, 
circular and elliptic geometry. Only one realization is chosen from each series. 

show ensemble averages of the absolute value of the normalized angular momentum over ten 
realizations. We observe that the magnitude of the spin-up increases more than a factor 2 
when increasing the Reynolds number by an order of magnitude. It is observed that the 
angular momentum in the circular domain is weaker but not negligible. 

In Figure |3] we show the angular momentum in the three geometries for series B and C 
corresponding to Reynolds numbers of order 10^ and 10^, respectively. For each Reynolds 
number one particular realization is chosen for which is maximum. For both series it is 
observed that strong spin-up takes place in the square and in the ellipse. The generation 
of the angular momentum is spontaneous, and rapid and one observes that the amplitude 
is of order 0.25 in the square and in the ellipse. This implies that the fluid reaches an 
angular momentum which corresponds to approximately 25% of the angular momentum 
which would possess a fluid in solid-body rotation containing the same energy at t = 0. 
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There is practically no spin-up in the circular container. 

In Figure 3, right, the magnetic angular momentum is evaluated in all geometries. Sur- 
prisingly, in the square in which the generation of kinetic angular momentum was the 
strongest. Lb remains close to zero. In the other two geometries, an amount of Lb is 
created, however, this magnetic spin-up takes place on a time-scale which is larger than 
for its kinetic counterpart. Furthermore it can be observed that once is created it re- 
mains almost constant over time. For series C Lb remains close to zero at all times in all 
geometries. 

C. Influence of the magnetic pressure 

IdIA -we derived the equation for L^ in the case of MHD turbulence. It reads 



with p the kinematic viscosity, u the vorticity, n the unit-vector perpendicular to the wall, 
p* = p + is the sum of the hydrodynamic and magnetic pressure. It was discovered by 
Clercx et al.^ that spontaneous generation of angular momentum in hydrodynamic turbu- 
lence is observed in square domains, whereas it is absent in a circular domain. Subsequently, 
it was explained to be an effect due to the pressure^, the last term in equation (fTT]) . Indeed, 
this term vanishes in a circular domain. In MHD, the presence of the magnetic pressure 
allows to vary the importance of the pressure term, while keeping the other parameters 
constant, by changing the value of the magnetic fluctuations. This is illustrated in Figure H] 
for series B (Reynolds ~ 10'^). The ratio Eb/E^ is varied, with Eb the mean-square of the 
magnetic fluctuations and E^ the mean-square of the velocity fluctuations. It is observed 
that the tendency to spin-up is significantly increased in the square geometry while this ef- 
fect is weaker in the elliptical geometry and absent in the circle. It is thus shown that both 
geometry and magnetic pressure can play a role in the generation of angular momentum. 

D. On the origin of the angular fields. 

In^'^, the tendency to generate angular fields was also investigated by computing the value 
of Lb- It was found that angular fields were observed, even in the circular geometry. In 
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Figure 3 right, we show that at higher Reynolds numbers the generation of this 'magnetic 
angular momentum' becomes weaker and seems to vanish. Writing the equation for L^, we 
find 

—f- = 7] (b j{r ■ n)ds - 2r//, (12) 
Jan 

where / denotes the net current through the domain, defined by J = J^jdA. The pressure 
plays thus no direct role and only the net current or resistive magnetic stress can generate 
angular fields. The mean current through the domain is computed by integrating the current 
density over the fluid domain. This quantity should in principle be small, and decay to zero 
at long times. No production of mean current is physically expected. Closer scrutiny of the 
results revealed the existence of a spurious fluctuating mean current inside the fluid domain. 
The fluctuations of this current are partly numerical. Indeed, the penalization method is 
known to induce small errors in the vicinity of the wall. These errors can be controlled and 
depend on the parameter e. The thickness of the layer in which the penalization error is 
significant is of order A = ^/eu. In this numerical boundary layer, non-physical currents can 
be observed. We will denote the total amount of numerical current by I^. If we suppose 
that this current is uniformly distributed in the boundary layer, we can write for a circular 
domain, 



2nrAjN (13) 



which gives an average numerical current density ~ InH'^t^'I'V^)- Now, equation ( 1T2 
becomes 



dL 



B 



dt 



rri2njN - 2rilN (14) 
-2v]in, (15) 



and for the special case of unity magnetic Prandtl number, v = r], this simplifies to 

dLs ( [v 



it - '''' 

The fact that we have a penalization parameter of the order of the viscosity leads to a 
non-negligible production of magnetic angular momentum through the dissipation term, 
proportional to 1^. As one can see in Figure [5l the time evolution of the mean current and 
the time derivative of the magnetic angular momentum, computed with a classical finite 
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difference scheme of first order, overlap quite well. Equation (fT6|) shows that the effect 
should become smaller when the ratio v/e is decreased. Since we used the same value for 
e in all runs and we decreased the viscosity to increase the Reynolds number, the influence 
of the current should become smaller at higher Reynolds number. Indeed in series C the 
generation of angular fields was dramatically reduced with respect to series B as observed 
in Figure |3l which confirms our assumption that the origin is due to a numerical boundary 
layer. A remaining open issue is why this effect was small or absent in the square geometry. 
We suspect that the effect is stronger for geometries in which the mask is not aligned with 
the numerical grid. Indeed, a so-called staircase-effect is expected to decrease the quality of 
the approximation near the walls. 

IV. CONCLUSIONS AND PERSPECTIVES 

In total 63 pseudo-spectral simulations of two-dimensional MHD turbulence in a bounded 
domain were performed. It was shown that spin-up takes place in non-axisymmetric geome- 
tries (squares, ellipses). This phenomenon, observed ir^ at low Reynolds number, persists 
at higher Reynolds numbers and becomes more pronounced. The generation of the magnetic 
equivalent of the angular momentum becomes much weaker at higher Reynolds numbers. 
The first effect, the kinetic spin- up can be enhanced by increasing the magnetic fluctuations. 
It is therefore clearly related to the pressure term p* . The generation of angular fields in 
our simulation was shown to have a numerical origin. The effect was argued to be related to 
the current density leaking into the domain and can therefore be physically relevant if the 
walls are not assumed to be insulated. Indeed, the influence of other boundary conditions 
constitute an interesting objective. The main objective remains however the investigation 
of the effect in fully three-dimensional unsteady MHD simulations. 
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FIG. 4: (Color online) Time evolution of angular momentum Lu{t) for series B {TZ ~ 10^). The 
influence of the magnetic pressure on the spin-up in the square, circle and ellipse is illustrated by 
changing the ratio Eb/E^, while keeping E^ fixed. The magnetic pressure is changed by varying 
Eb, while keeping constant E^- 
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FIG. 5: (Color online) Comparison of the time derivative of Lsit) and the mean current < I >. 
The run corresponds to one realization in the circle with Eb/Eu = 13.3 and TZ ~ 10^. 
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